#clean up
rm(list=ls())

#load libraries
library(lattice)
library(xtable)
library(foreign)
library(maps)
library(geoR)
library(rgdal)
library(car)
library(maptools)

#options
options(scipen=8)

#set to YOUR working directory
#setwd('/Users/jamie/looseData/krigingOutput/')

#### PROCESSING MCMC SAMPLES ####
for(i in 1:100){
	check<-file.exists(paste('sample',i,'.csv',sep=''))
	if(check){
		assign(paste('s',i,sep=''),read.csv(paste('sample',i,'.csv',sep='')))
	}else{
		print(paste('File sample',i,'.csv does not exist.',sep=''))
		}
	}

#merge and drop component files
mcmc<-rbind(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,s15,s16,s17,s18,s19,s20,s21,s22,s23,s24,s25,s26,s27,s28,s29,s30,s31,s32,s33,s34,s35,s36,s37,s38,s39,s40,s41,s42,s43,s44,s45,s46,s47,s48,s49,s50,s51,s52,s53,s54,s55,s56,s57,s58,s59,s60,s61,s62,s63,s64,s65,s66,s67,s68,s69,s70,s71,s72,s73,s74,s75,s76,s77,s78,s79,s80,s81,s82,s83,s84,s85,s86,s87,s88,s89,s90,s91,s92,s93,s94,s95,s96,s97,s98,s99,s100)
rm(s1,s2,s3,s4,s5,s6,s7,s8,s9,s10,s11,s12,s13,s14,s15,s16,s17,s18,s19,s20,s21,s22,s23,s24,s25,s26,s27,s28,s29,s30,s31,s32,s33,s34,s35,s36,s37,s38,s39,s40,s41,s42,s43,s44,s45,s46,s47,s48,s49,s50,s51,s52,s53,s54,s55,s56,s57,s58,s59,s60,s61,s62,s63,s64,s65,s66,s67,s68,s69,s70,s71,s72,s73,s74,s75,s76,s77,s78,s79,s80,s81,s82,s83,s84,s85,s86,s87,s88,s89,s90,s91,s92,s93,s94,s95,s96,s97,s98,s99,s100)
#tau squared alone
mcmc$tausq.alone<-mcmc$tausq.rel*mcmc$sigmasq

#Density Plots
pdf('tausq.pdf',width=5,height=5);densityplot(mcmc$tausq.alone,xlab=expression(tau^2));dev.off()
pdf('tausqRel.pdf',width=5,height=5);densityplot(mcmc$tausq.rel,xlab=expression(tau^2/sigma^2));dev.off()
pdf('phi.pdf',width=5,height=5);densityplot(mcmc$phi,xlab=expression(1/phi));dev.off()
pdf('sigmasq.pdf',width=5,height=5);densityplot(mcmc$sigmasq,xlab=expression(sigma^2));dev.off()
postscript('tausq.eps',width=5,height=5);densityplot(mcmc$tausq.alone,xlab=expression(tau^2));dev.off()
postscript('tausqRel.eps',width=5,height=5);densityplot(mcmc$tausq.rel,xlab=expression(tau^2/sigma^2));dev.off()
postscript('phi.eps',width=5,height=5);densityplot(mcmc$phi,xlab=expression(1/phi));dev.off()
postscript('sigmasq.eps',width=5,height=5);densityplot(mcmc$sigmasq,xlab=expression(sigma^2));dev.off()

#rescale eastings and northings polynomial coefficients
mcmc$beta30<-mcmc$beta30*1000
mcmc$beta31<-mcmc$beta31*1000
mcmc$beta32<-mcmc$beta32*1000000
mcmc$beta33<-mcmc$beta33*1000000
mcmc$beta37<-mcmc$beta37*1000000

#Output in a LaTeX table
xtable(
cbind(
apply(mcmc,2,mean),
#apply(mcmc,2,median),
apply(mcmc,2,sd),
apply(mcmc,2,quantile,.05),
apply(mcmc,2,quantile,.95)
),
digits=4)

###load original data###
combined<-read.dta('cces08reformatted.dta')

###plot empirical variogram###
#create variograms
ideol.robust<-variog(coords=cbind(combined$eastings,combined$northings), data=combined$ideology, estimator.type="modulus")
ideol.ols<-lm(ideology~age*educ+as.factor(race)*female+as.factor(inc14)+catholic+mormon+orthodox+jewish+islam+mainline+evangelical+rural+ownership+as.factor(empstat)+eastings*northings+I(eastings^2)+I(northings^2),data=combined);summary(ideol.ols)
combined$ols.resid<-ideol.ols$residuals
resid.robust<-variog(coords=cbind(combined$eastings,combined$northings), data=combined$ols.resid, estimator.type="modulus")
gaussian.fit <- variofit(resid.robust, ini.cov.pars = c(100,3000), cov.model="gaussian", fix.nugget=FALSE, nugget=700)

#OLS output
#postscript('olsGaussian.eps',width=5,height=5)
#pdf('olsGaussian.pdf',width=5,height=5)
plot(ideol.robust, xlab="Distance in Kilometers", ylab="Semivariance",ylim=c(600,725))
par(new=T)
plot(resid.robust, xlab="", ylab="",axes=F,ylim=c(600,800),pch='+', col='blue')
lines(gaussian.fit, col='blue')
#abline(v=384144)
#dev.off()

#MCMC Output
#postscript('mcmcGaussian.eps',width=5,height=5)
#pdf('mcmcGaussian.pdf',width=5,height=5)
plot(ideol.robust, xlab="Distance in Kilometers", ylab="Semivariance",ylim=c(600,725))
par(new=T)
plot(resid.robust, xlab="", ylab="",axes=F,ylim=c(600,800),pch='+', col='blue')
lines.variomodel(cov.model="gaussian",cov.pars=c(86.8419,5007.6302),nug=604.1449,col='red',max.dist=5000)
#abline(v=3602.04)
#dev.off()

#### PROCESSING FORECASTS ####
for(i in 1:100){
	check<-file.exists(paste('forecast',i,'.csv',sep=''))
	if(check){
		assign(paste('f',i,sep=''),read.csv(paste('forecast',i,'.csv',sep='')))
	}else{
		print(paste('File forecast',i,'.csv does not exist.',sep=''))
		}
	}

#merge and drop component files
forecasts<-rbind(f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18,f19,f20,f21,f22,f23,f24,f25,f26,f27,f28,f29,f30,f31,f32,f33,f34,f35,f36,f37,f38,f39,f40,f41,f42,f43,f44,f45,f46,f47,f48,f49,f50,f51,f52,f53,f54,f55,f56,f57,f58,f59,f60,f61,f62,f63,f64,f65,f66,f67,f68,f69,f70,f71,f72,f73,f74,f75,f76,f77,f78,f79,f80,f81,f82,f83,f84,f85,f86,f87,f88,f89,f90,f91,f92,f93,f94,f95,f96,f97,f98,f99,f100)

rm(f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11,f12,f13,f14,f15,f16,f17,f18,f19,f20,f21,f22,f23,f24,f25,f26,f27,f28,f29,f30,f31,f32,f33,f34,f35,f36,f37,f38,f39,f40,f41,f42,f43,f44,f45,f46,f47,f48,f49,f50,f51,f52,f53,f54,f55,f56,f57,f58,f59,f60,f61,f62,f63,f64,f65,f66,f67,f68,f69,f70,f71,f72,f73,f74,f75,f76,f77,f78,f79,f80,f81,f82,f83,f84,f85,f86,f87,f88,f89,f90,f91,f92,f93,f94,f95,f96,f97,f98,f99,f100)

#obtain state averages from kriging
stateKriging<-aggregate(x=forecasts$predictive.mean, by=list(FIPS=forecasts$statefip), FUN=mean)
colnames(stateKriging)[2]<-'krige.state'
between<-aggregate(x=forecasts$predictive.mean, by=list(FIPS=forecasts$statefip), FUN=var)
within<-aggregate(x=forecasts$predictive.variance, by=list(FIPS=forecasts$statefip), FUN=mean)
no.kriges<-table(forecasts$statefip)
stateKriging$krige.state.var<-within[,2]+((no.kriges+1)/no.kriges)*between[,2]

#write output
write.csv(stateKriging, 'stateKriging.csv', row.names=F)

###View the performance of the forecasts
#NEW working directory. Set to YOUR working directory.
#setwd('/Volumes/MONOGAN/gillAreal/ccesCode/')

#load data from forecasts and elections
stateMeasures<-read.dta('stateMeasures.dta')
cor(stateMeasures$prespct, stateMeasures$krigestate)	#r=-0.8280817
cor(stateMeasures$prespct, stateMeasures$krigestate)^2	#r^2=0.6857194
cor(stateMeasures$prespct, stateMeasures$unweighted) #r=-0.8911825
cor(stateMeasures$prespct, stateMeasures$weighted) #r=-0.8344995
cor(stateMeasures$prespct, stateMeasures$mrpstate) #r=-0.8667577

postscript('statePres.eps',width=5,height=5)
#pdf('statePres.pdf',width=5,height=5)
plot(y=stateMeasures$prespct, x=stateMeasures$krigestate, type='n',xlab='Kriging Forecast of Ideology',ylab="Obama's Vote Share (2008)")
text(y=stateMeasures$prespct, x=stateMeasures$krigestate, labels=as.character(stateMeasures$state))
abline(reg=lm(stateMeasures$prespct~stateMeasures$krigestate),col='red')
dev.off()

postscript('stateUnweighted.eps',width=5,height=5)
#pdf('stateSubsetting.pdf',width=5,height=5)
plot(y=stateMeasures$prespct, x=stateMeasures$unweighted, type='n',xlab='State Subsets of Ideology',ylab="Obama's Vote Share (2008)")
text(y=stateMeasures$prespct, x=stateMeasures$unweighted, labels=as.character(stateMeasures$state))
abline(reg=lm(stateMeasures$prespct~stateMeasures$unweighted),col='red')
dev.off()

postscript('stateMRP.eps',width=5,height=5)
#pdf('stateMRP.pdf',width=5,height=5)
plot(y=stateMeasures$prespct, x=stateMeasures$mrpstate, type='n',xlab='MRP Forecast of Ideology',ylab="Obama's Vote Share (2008)")
text(y=stateMeasures$prespct, x=stateMeasures$mrpstate, labels=as.character(stateMeasures$state))
abline(reg=lm(stateMeasures$prespct~stateMeasures$mrpstate),col='red')
dev.off()

#make a map of ideology, as we measured it
X <- c(stateMeasures$krigestate[1:20],stateMeasures$krigestate[20],stateMeasures$krigestate[20:21],stateMeasures$krigestate[21:31],stateMeasures$krigestate[31],stateMeasures$krigestate[31],stateMeasures$krigestate[31:32],stateMeasures$krigestate[32],stateMeasures$krigestate[32:45],stateMeasures$krigestate[45],stateMeasures$krigestate[45:46],stateMeasures$krigestate[46],stateMeasures$krigestate[46],stateMeasures$krigestate[46],stateMeasures$krigestate[46:49])
n.col <- 10
br <- c(-999,quantile(X,c(.1,.2,.3,.4,.5,.6,.7,.8,.9)), 999)
#shading<-gray((n.col-1):0/(n.col-1))
shading<- colorRampPalette(c("blue","white","red"))(10)
#shading<- colorRampPalette(c("red","blue"))(10)
X.grp<-findInterval(X, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.shad<-shading[X.grp]
IDs <- map("state",fill=FALSE, plot=FALSE, resolution=0)$names

postscript('krigingMap.eps',width=5,height=5)
#pdf('krigingMap.pdf',width=5,height=5)
map(database="state", region=IDs, fill=T, plot=T, interior=T, col=X.shad, exact=T, resolution=0)
dev.off()

#make a map of ideology, using subsets
X <- c(stateMeasures$unweighted[1:20],stateMeasures$unweighted[20],stateMeasures$unweighted[20:21],stateMeasures$unweighted[21:31],stateMeasures$unweighted[31],stateMeasures$unweighted[31],stateMeasures$unweighted[31:32],stateMeasures$unweighted[32],stateMeasures$unweighted[32:45],stateMeasures$unweighted[45],stateMeasures$unweighted[45:46],stateMeasures$unweighted[46],stateMeasures$unweighted[46],stateMeasures$unweighted[46],stateMeasures$unweighted[46:49])
n.col <- 10
br <- c(-999,quantile(X,c(.1,.2,.3,.4,.5,.6,.7,.8,.9)), 999)
#shading<-gray((n.col-1):0/(n.col-1))
shading<- colorRampPalette(c("blue","white","red"))(10)
#shading<- colorRampPalette(c("red","blue"))(10)
X.grp<-findInterval(X, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.shad<-shading[X.grp]
IDs <- map("state",fill=FALSE, plot=FALSE, resolution=0)$names

postscript('subsetMap.eps',width=5,height=5)
#pdf('subsetMap.pdf',width=5,height=5)
map(database="state", region=IDs, fill=T, plot=T, interior=T, col=X.shad, exact=T, resolution=0)
dev.off()

#make a map of MRP state ideology
X <- c(stateMeasures$mrpstate[1:20],stateMeasures$mrpstate[20],stateMeasures$mrpstate[20:21],stateMeasures$mrpstate[21:31],stateMeasures$mrpstate[31],stateMeasures$mrpstate[31],stateMeasures$mrpstate[31:32],stateMeasures$mrpstate[32],stateMeasures$mrpstate[32:45],stateMeasures$mrpstate[45],stateMeasures$mrpstate[45:46],stateMeasures$mrpstate[46],stateMeasures$mrpstate[46],stateMeasures$mrpstate[46],stateMeasures$mrpstate[46:49])
n.col <- 10
br <- c(-999,quantile(X,c(.1,.2,.3,.4,.5,.6,.7,.8,.9)), 999)
#shading<-gray((n.col-1):0/(n.col-1))
shading<- colorRampPalette(c("blue","white","red"))(10)
#shading<- colorRampPalette(c("red","blue"))(10)
X.grp<-findInterval(X, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.shad<-shading[X.grp]
IDs <- map("state",fill=FALSE, plot=FALSE, resolution=0)$names

postscript('mrpMap.eps',width=5,height=5)
#pdf('mrpMap.pdf',width=5,height=5)
map(database="state", region=IDs, fill=T, plot=T, interior=T, col=X.shad, exact=T, resolution=0)
dev.off()

#make a map of kriging variance
X <- c(stateMeasures$krigestatevar[1:20],stateMeasures$krigestatevar[20],stateMeasures$krigestatevar[20:21],stateMeasures$krigestatevar[21:31],stateMeasures$krigestatevar[31],stateMeasures$krigestatevar[31],stateMeasures$krigestatevar[31:32],stateMeasures$krigestatevar[32],stateMeasures$krigestatevar[32:45],stateMeasures$krigestatevar[45],stateMeasures$krigestatevar[45:46],stateMeasures$krigestatevar[46],stateMeasures$krigestatevar[46],stateMeasures$krigestatevar[46],stateMeasures$krigestatevar[46:49])
n.col <- 10
br <- c(-999,quantile(X,c(.1,.2,.3,.4,.5,.6,.7,.8,.9)), 999)
#shading<-gray((n.col-1):0/(n.col-1))
shading<- colorRampPalette(c("green","white","orange"))(10)
#shading<- colorRampPalette(c("red","blue"))(10)
X.grp<-findInterval(X, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.shad<-shading[X.grp]
IDs <- map("state",fill=FALSE, plot=FALSE, resolution=0)$names

postscript('krigingVariance.eps',width=5,height=5)
#pdf('krigingVariance.pdf',width=5,height=5)
map(database="state", region=IDs, fill=T, plot=T, interior=T, col=X.shad, exact=T, resolution=0)
dev.off()

#make a map of presidential vote
X <- c(stateMeasures$prespct[1:20],stateMeasures$prespct[20],stateMeasures$prespct[20:21],stateMeasures$prespct[21:31],stateMeasures$prespct[31],stateMeasures$prespct[31],stateMeasures$prespct[31:32],stateMeasures$prespct[32],stateMeasures$prespct[32:45],stateMeasures$prespct[45],stateMeasures$prespct[45:46],stateMeasures$prespct[46],stateMeasures$prespct[46],stateMeasures$prespct[46],stateMeasures$prespct[46:49])
n.col <- 10
br <- c(-999,quantile(X,c(.1,.2,.3,.4,.5,.6,.7,.8,.9)), 999)
#shading<-gray((n.col-1):0/(n.col-1))
shading<- colorRampPalette(c("red","white","blue"))(10)
#shading<- colorRampPalette(c("red","blue"))(10)
X.grp<-findInterval(X, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.shad<-shading[X.grp]
IDs <- map("state",fill=FALSE, plot=FALSE, resolution=0)$names

postscript('statePresMap.eps',width=5,height=5)
#pdf('statePresMap.pdf',width=5,height=5)
map(database="state", region=IDs, fill=T, plot=T, interior=T, col=X.shad, exact=T, resolution=0)
dev.off()



#### PROCESSING DISTRICT-LEVEL FORECASTS ####
#merge and drop component files
vaDist<-subset(forecasts,!is.na(cd))

#obtain state averages from kriging
distKriging<-aggregate(x=vaDist$predictive.mean, by=list(cd=vaDist$cd), FUN=mean)
colnames(distKriging)[2]<-'krige.cd'
between<-aggregate(x=vaDist$predictive.mean, by=list(cd=vaDist$cd), FUN=var)
within<-aggregate(x=vaDist$predictive.variance, by=list(cd=vaDist$cd), FUN=mean)
no.kriges<-table(vaDist$cd)
distKriging$krige.cd.var<-within[,2]+((no.kriges+1)/no.kriges)*between[,2]
write.csv(distKriging, 'distKriging.csv', row.names=F)

#load map of VA congressional districts
va.cd.0<-readShapePoly(fn="/Volumes/MONOGAN/gillAreal/ccesCode/tl_2008_51_cd110/tl_2008_51_cd110.shp", proj4string=CRS("+proj=longlat +datum=NAD83"))
va.cd.1<-spTransform(va.cd.0,CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=km +no_defs"))
va.cd.2<-as(va.cd.1,'SpatialPolygons')

###View the performance of the forecasts
distMeasures<-read.dta('vaDistMeasures.dta')
cor(distMeasures$obama08, distMeasures$krigecd)	#r=-0.774952
cor(distMeasures$obama08, distMeasures$krigecd)^2	#r^2=0.6005506
cor(distMeasures$dem08[-c(3,9)], distMeasures$krigecd[-c(3,9)])	#r=-0.6269749

cor(distMeasures$obama08, distMeasures$unweighted)	#r=-0.6001
cor(distMeasures$obama08, distMeasures$unweighted)^2	#r^2=0.3601
cor(distMeasures$dem08[-c(3,9)], distMeasures$unweighted[-c(3,9)])	#r=-0.3213

cor(distMeasures$obama08, distMeasures$weighted)	#r=-0.4934
cor(distMeasures$obama08, distMeasures$weighted)^2	#r^2=0.2435
cor(distMeasures$dem08[-c(3,9)], distMeasures$weighted[-c(3,9)])	#r=-0.2823

cor(distMeasures$obama08, distMeasures$mrpcd)	#r=-0.8852459
cor(distMeasures$obama08, distMeasures$mrpcd)^2	#r^2=0.7836603
cor(distMeasures$dem08[-c(3,9)], distMeasures$mrpcd[-c(3,9)])	#r=-0.6017974

postscript('distPres.eps',width=5,height=5)
#pdf('distPres.pdf',width=5,height=5)
plot(y=distMeasures$obama08, x=distMeasures$krigecd, pch=19,xlab='Kriging Forecast of Ideology',ylab="Obama's Vote Share (2008)")
abline(reg=lm(distMeasures$obama08~distMeasures$krigecd),col='red')
dev.off()

postscript('distSubsetting.eps',width=5,height=5)
#pdf('distSubsetting.pdf',width=5,height=5)
plot(y=distMeasures$obama08, x=distMeasures$unweighted, pch=19,xlab='District Subsets of Ideology',ylab="Obama's Vote Share (2008)")
abline(reg=lm(distMeasures$obama08~distMeasures$unweighted),col='red')
dev.off()

postscript('distMRP.eps',width=5,height=5)
#pdf('distMRP.pdf',width=5,height=5)
plot(y=distMeasures$obama08, x=distMeasures$mrpcd, pch=19,xlab='MRP Forecast of Ideology',ylab="Obama's Vote Share (2008)")
abline(reg=lm(distMeasures$obama08~distMeasures$mrpcd),col='red')
dev.off()

postscript('congKrig.eps',width=5,height=5)
#pdf('distPres.pdf',width=5,height=5)
plot(y=distMeasures$dem08[-c(3,9)], x=distMeasures$krigecd[-c(3,9)], pch=19,xlab='Kriging Forecast of Ideology',ylab="Democratic Vote Share, House Races (2008)")
abline(reg=lm(distMeasures$dem08[-c(3,9)]~distMeasures$krigecd[-c(3,9)]),col='red')
dev.off()

postscript('congSubset.eps',width=5,height=5)
#pdf('distSubsetting.pdf',width=5,height=5)
plot(y=distMeasures$dem08[-c(3,9)], x=distMeasures$unweighted[-c(3,9)], pch=19,xlab='District Subsets of Ideology',ylab="Democratic Vote Share, House Races (2008)")
abline(reg=lm(distMeasures$dem08[-c(3,9)]~distMeasures$unweighted[-c(3,9)]),col='red')
dev.off()

postscript('congMRP.eps',width=5,height=5)
#pdf('congMRP.pdf',width=5,height=5)
plot(y=distMeasures$dem08, x=distMeasures$mrpcd, pch=19,xlab='MRP Forecast of Ideology',ylab="Obama's Vote Share (2008)")
abline(reg=lm(distMeasures$dem08~distMeasures$mrpcd),col='red')
dev.off()


#make a map of district ideology, as we measured it
X <- (distMeasures$krigecd)
n.col <- 10
br <- c(-999,quantile(X,c(.1,.2,.3,.4,.5,.6,.7,.8,.9)), 999)
#shading<-gray((n.col-1):0/(n.col-1))
shading<- colorRampPalette(c("blue","white","red"))(10)
#shading<- colorRampPalette(c("red","blue"))(10)
X.grp<-findInterval(X, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.shad<-shading[X.grp][c(6,9,10,1,11,5,7,8,4,3,2)]

#Charlottesville
#charlottesville.0<-SpatialPoints(coords=matrix(c(-78.4865, 38.0345),nrow=1),proj4string=CRS("+proj=longlat +datum=NAD83"))
#charlottesville<-spTransform(charlottesville.0,CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=km +no_defs"))

postscript('vaKrigeMap.eps',width=5,height=5)
#pdf('vaKrigeMap.pdf',width=5,height=5)
plot(va.cd.2, col=X.shad)
#points(charlottesville,pch=15)
dev.off()

#make a map of presidential vote
X <- (distMeasures$obama08)
n.col <- 10
br <- c(-999,quantile(X,c(.1,.2,.3,.4,.5,.6,.7,.8,.9)), 999)
#shading<-gray((n.col-1):0/(n.col-1))
shading<- colorRampPalette(c("red","white","blue"))(10)
#shading<- colorRampPalette(c("red","blue"))(10)
X.grp<-findInterval(X, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.shad<-shading[X.grp][c(6,9,10,1,11,5,7,8,4,3,2)]

postscript('vaPresMap.eps',width=5,height=5)
#pdf('vaPresMap.pdf',width=5,height=5)
plot(va.cd.2, col=X.shad)
#points(charlottesville,pch=15)
dev.off()

#make a map of district ideology, using subsets
X <- (distMeasures$unweighted)
n.col <- 10
br <- c(-999,quantile(X,c(.1,.2,.3,.4,.5,.6,.7,.8,.9)), 999)
#shading<-gray((n.col-1):0/(n.col-1))
shading<- colorRampPalette(c("blue","white","red"))(10)
#shading<- colorRampPalette(c("red","blue"))(10)
X.grp<-findInterval(X, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.shad<-shading[X.grp][c(6,9,10,1,11,5,7,8,4,3,2)]

postscript('vaSubsetMap.eps',width=5,height=5)
#pdf('vaSubsetMap.pdf',width=5,height=5)
plot(va.cd.2, col=X.shad)
#points(charlottesville,pch=15)
dev.off()

#make a map of district ideology, using MRP
X <- (distMeasures$mrpcd)
n.col <- 10
br <- c(-999,quantile(X,c(.1,.2,.3,.4,.5,.6,.7,.8,.9)), 999)
#shading<-gray((n.col-1):0/(n.col-1))
shading<- colorRampPalette(c("blue","white","red"))(10)
#shading<- colorRampPalette(c("red","blue"))(10)
X.grp<-findInterval(X, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.shad<-shading[X.grp][c(6,9,10,1,11,5,7,8,4,3,2)]

postscript('vaMRPmap.eps',width=5,height=5)
#pdf('vaMRPmap.pdf',width=5,height=5)
plot(va.cd.2, col=X.shad)
#points(charlottesville,pch=15)
dev.off()
